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Abstract. While there are many well-known and extensively tested results involving 
diffusion-limited binary reactions, reactions involving subdiffusive reactant species 
are far less understood. Subdiffusive motion is characterized by a mean square 
displacement (x^) ~ f with < 7 < 1. Recently we calculated the asymptotic survival 
probability P{t) of a (sub)diffusive particle (7') surrounded by (sub)diffusive traps (7) 
in one dimension. These are among the few known results for reactions involving 
species characterized by different anomalous exponents. Our results were obtained by 
bounding, above and below, the exact survival probability by two other probabilities 
that are asymptotically identical (except when 7' = 1 and < 7 < 2/3). Using this 
approach, we were not able to estimate the time of validity of the asymptotic result, 
nor the way in which the survival probability approaches this regime. Toward this 
goal, here we present a detailed comparison of the asymptotic results with numerical 
simulations. In some parameter ranges the asymptotic theory describes the simulation 
results very well even for relatively short times. However, in other regimes more time is 
required for the simulation results to approach asymptotic behavior, and we arrive at 
situations where we are not able to reach asymptotia within our computational means. 
This is regrettably the case for 7' = 1 and < 7 < 2/3, where we are therefore not 
able to prove or disprove even conjectures about the asymptotic survival probability 
of the particle. 
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1. Introduction 

The survival probability of a particle diffusing in a one-dimensional medium of diffusive 
traps has only recently been calculated (but only asymptotically) [1-5]. This is 
surprising in view of its long history [6-8] and that of its antecedents, the so-called 
trapping problem [9-18], in which the traps are static, and the target problem [15-18], 
in which the particle does not move. The antecedent systems could be translated to 
tractable boundary value problems, which is not possible when both particle and traps 
move. The solution is an elegant "tour de force" in which the desired survival probability 
is bounded above and below by two others that can be posed as boundary value problems 
and that converge to one another asymptotically. 

Recently, we undertook the generalization of the bounding approach to the case 
of a subdiffusive particle surrounded by a distribution of subdiffusive traps [19,20]. 
Subdiffusion of a particle is usually characterized by the time dependence of the mean 
square of the particle displacement x(t), 

Here is the (generalized) diffusion constant, and 7 is the exponent that characterizes 
normal (7 = 1) or anomalous (7 7^ 1) diffusion. In particular, the process is diffusive 
when 7 = 1 and sudiffusive when < 7 < 1. There are a variety of models and physical 
circumstances that lead to subdiffusion in the trapping and, more generally, in the 
binary reaction context [15-17,21-31]. Many are based on the continuous time random 
walk formalism, where particles are thought of as random walkers with waiting-time 
distributions between steps that have broad long-time tails and consequently infinite 
moments, ip{t) ~ t~^~'^ . Our work is based on the fractional diffusion equation, which 
describes the evolution of the probability density P{x, t) of finding the particle at 
position X at time t by means of the fractional partial differential equation (in one 
dimension) [32,33], 

^^P{x,t) = K,oDl-'^P{x,t), (2) 

1-7 



where is the generalized diffusion coefficient that appears in equation ([T]), and -D- 
is the Riemann-Liouville operator. 



oDr'P{x,t) = -—^ dr ^ ' . (3) 

The connection between these two approaches is in itself an interesting subject, see 
e.g. [34,35]. 

The survival probability of a particle A characterized by exponent 7' and generalized 
diffusion coefficient Ky surrounded by traps characterized by 7 and is bounded as 
follows [19,20]. An upper bound is obtained by forcing particle A to remain still. The 
"Pascal principle" that says that the best survival strategy for the particle is to stand 
still was proved for the diffusive case in [4, 5, 8] and for the subdiffusive problem in 
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our work. The solution of the fractional subdiffusion equation for the B traps with the 
location of A as an appropriate boundary then leads to the upper bound for the survival 
probability of A, 



Pu{t) = exp 



-2p 



K^t-r 



r(i + 



(4) 



where p is the density of traps. A lower bound is calculated by allowing particle A to 
move within a box of size C while the traps B are forced to remain outside of this box. 
The box size is then found so as to maximize this lower bound, with the t oo result 

2 




X exp 



1 + 



o'K'ty 



(5) 



for < 7' < 1. For the diffusive case (7' 

-2p 



4 



PL{t) = - exp 



TT 



K^t-y 



r(i + 7/2) 



3(7rVD't/4) 



1/3 



(6) 



With these bounds we see that for a subdiffusive particle (0 < 7' < 1) and diffusive 
or subdiffusive traps (0 < 7 < 1) the upper and lower bounds converge asymptotically 



(viz., compare the logarithms of both), 
survival probability 



P{t) ~ exp 



-2p 



r(i 



SO that we arrive at the explicit asymptotic 



(7) 



This result elicits a comment about the so-called "subordination principle" [17], 
according to which in some cases asymptotic anomalous diffusion behavior can be found 
from corresponding results for normal diffusion with the simple replacement of t by 
fy . This can be understood from a continuous time random walk perpective because 
the average number of jumps made by a subdiffusive walker up to time t scales as 
(n) ~ fy, and in many instances the number of jumps is the relevant factor that 
explains the behavior of the system. However, for systems where each species has a 
different anomalous diffusion exponent, such a replacement becomes ambiguous. The 
result ([7]) indicates a subordination principle at work as determined by the traps. In 
other words, it is the motion of the traps that regulates the survival probability of the 
particle whether or not the particle moves, provided it does not move "too easily," i.e., 
provided it is subdiffusive (7' < 1). 

When the particle is diffusive (7' = 1) the situation is more complicated, because 
its asymptotic survival probability is no longer necessarily the same as it would be if it 
stood still. If 2/3 < 7 < 1, i.e., if the traps move sufficiently easily, the upper and lower 
bounds still converge and Eq. ([7]) still holds, that is, it is still the motion of the traps 
that determines the asymptotic survival probability of the particle. If the traps are 
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sub diffusive with 7 = 2/3 (marginal case), the bounds lead only to a prediction of the 
asymptotic time dependence but not of the accompanying exponential prefactor, i.e., the 
bounds establish that P{t) ~ exp(— At^/^) but are not able to determine A. In particular. 



we can not determine whether whether A is given by the coefficient 2py^K2/3/r(4/3) of 

t^/^ in the exponent of Eq. ([7]) when 7 = 2/3, which one might conjecture. Finally, when 
the particle is diffusive and the traps are sufficiently slow (0 < 7 < 2/3), the upper and 
lower bounds do not have the same asymptotic time dependence, so we are not able to 
even functionally bound the survival probability. While it is still possible in principle 
that the trap-driven subordination principle continues to apply in this regime so that 
6 = 7/2, we have not been able to prove or disprove such a conjecture (and the behavior 
at 7 = would not fall within this conjecture, see below). If this subordination result 
is invalid, it would imply (and would not be surprising) that it is no longer possible to 
assume the particle to be standing still, and/or that it may no longer be (or only be) the 
exponent of the traps that regulates the motion. It is interesting to note that for 7' = 1 
and 7 = 0, the traditional "trapping problem," the asymptotic survival probability 
is [36] P{t) ~ exp[— 3(7r^p^D/4)^/^t^/^]. One might thus be tempted to conjecture a 
behavior of the form P(t) ~ exp(— A(7, p)t^/^) throughout the range < 7 < 2/3, i.e., 
the exponential prefactor A would have to depend on 7 and on p. However, we have not 
been able to prove or disprove this conjecture either. 

This analysis therefore leaves open two important questions, which we attempt to 
answer by way of detailed numerical simulations (although we do not entirely succeed): 

(i) In the cases where the bounds converge asymptotically, how much time does it take 
for the result to adequately describe the survival probability of the particle? In 
other words, how rapidly do the upper and lower bounds converge to one another? 

(ii) Is it possible to find the asymptotic survival probability for the cases in which our 
analysis fails to provide converging bounds? 

In part the success or failure of this attempt is of course constrained by the numerical 
resources at our disposal; the difficulties in numerically reaching asymptotia even in 
diffusion problems are well known [37]. 

The rest of this paper thus consists mostly of figures and a table presenting 
numerical simulation results. Our purpose is to ascertain the behavior of A and 6 
in the expression 



if indeed we arrive at a regime where the survival probability exhibits this behavior. 
If asymptotia has been reached, we expect both to be constant with time. If they are 
constant, we compare their values with those obtained from the asymptotic result ([7]). 
To test 6, we plot ln[— InP(t)] vs Int. For some cases where 6 is clearly determined as 
a result of our simulations (i.e., constant in time and independent of p), we test A by 
plotting -\nP{t)/x vs x, 




P(t) = exp(-Ar) 



(8) 




(9) 
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X is in effect a convenient dimensionless measure of time, and is also the ratio of tlie root 
mean square displacement of a particle at time t to the average distance between 
traps. 

In section |2] we briefly outline our simulation methodology. Section [3] is a 
compendium of our results, along with the associated descriptions. A recapitulation 
is presented in section HJ 

2. Numerical simulation methodology 

A brief review of our numerical simulation methodology is appropriate at this point. We 
generate the trap distribution by placing a trap at each site of a one-dimensional lattice 
with probability p (and not placing a trap with probability 1 — p). The particle A is 
placed at the origin of the lattice. The typical lattice has 10 000 sites, and we implement 
periodic boundary conditions. We have simulated larger lattices and different (free) 
boundary conditions to ascertain that the results are not affected. 

The dynamics of a moving particle in a sea of moving traps is implemented as 
follows. Each particle and trap is assigned an "internal clock" starting at time t = 
according to their waiting time probability distributions. One particular trap, or the 
particle, will be the first to take a step, left or right with equal probability (1/2). We 
check if trapping of the particle occurs as a result. If it does, we stop the dynamics, 
record the time, and generate a new ensemble of traps plus one particle. If it does 
not, we continue the dynamics by observing the very next trap or particle that takes 
a step. Again, if trapping occurs, the time is recorded and the dynamics stopped; if 
not, the walk continues. We also define a maximal time threshold (dictated by our 
computational resources) at which we stop the dynamics. 

In order to collect enough statistics we have run a large number of realizations 
(ensembles of traps and particle) of the dynamics, typically on the order of 50 000. A 
64-bit congruential random number generator was used througout the program [38]. The 
output of interest of each realization is the time when the particle is annhilated. On the 
basis of this observable we construct the integrated probability distribution of particle 
survival. The statistical errors have been computed using the jackknife procedure [39]. 

3. Results 

We will see that our numerical simulations of the problem A + B ^ B when both the 
traps B and the particle A are mobile approach asymptotic behavior far more rapidly in 
some regions of parameter space than in others. Figured] shows a 7-7' space in which are 
indicated regions and points of exponents 6 predicted analytically. The thick (orange 
in colour rendition) strip on the right (not including its end points) represents the 
parameter regime where the upper and lower survival probability bounds do not coverge 
asymptotically and hence no predictions (aside from conjectures) have been made. The 
prefactor A has also been predicted everywhere except on the thick (orange) strip and 
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Figure 1. 7-7' parameter space indicating regions of analytic asymptotic predictions 
for the exponent 9 and the coefficient A in P{t) ^ exp(— At^). No predictions are 
available along the thick strip (orange in colour version) 7' = 1, < 7 < 2/3. The 
exponent is predicted at the end points of the thick strip, the prefactor only at the 
lower end point. 



its upper end point, indicated by a circle (green in colour rendition). The empty squares 
in figure [1] indicate the pairs (7', 7) where we have carried out numerical simulations. 
All of our results for the apparent exponent 6 obtained from the simulations for these 
points are summarized in table [H and a number of them are subsequently exhibited in 
figures. 

The table leads to a number of broad conclusions, starting with the assertion that 
in the parameter regime 7 > 7' we have been able to reach the asymptotic exponent 6, 
and that this exponent agrees with the theoretical asymptotic prediction. In this regime 
the slope we identify as 9 is indeed insensitive to trap density changes and close to the 
asymptotic value 7/2. We have listed four sets of results for 7 > 7', and exhibit three 
of them explicitly in figures [21 [3], and |H We see that the time to achieve asymptotic 
behavior increases drastically with decreasing 7 for a given p. For example, for the same 
p in figures [2] and [3], we need to go to ln(t) > 20 or so for 7 = 0.4 but only to ln(t) > 10 
or so for 7 = 0.8. For a given pair of parameters (7, 7'), the time needed to reach 
asymptotic behavior is of course greater when the density of traps is lower. Note also 
that this time seems insensitive to the value of 7': for the case in the table not shown 
in a figure (7 = 7' = 0.8), the abscissa ln(t) covers the same range for the same trap 
densities as those in figure [31 

The situation is more complicated and less satisfactory for 7 < 7'. This is seen into 
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Figure 2. Simulation results for 7 = 7' = 0.4. The left set of data (black in colour 
version) corresponds to a trap density p = 0.5 and a slope of 0.210; the right set (red 
in the colour version) to p = 0.1 and a slope of 0.186. The asymptotic prediction for 
the slope is 7/2 = 0.2. The error bars of the points are smaller than the size of the 
symbols (the same in figures [3] to [8]). 




ln(t) 



Figure 3. Simulation results for 7 = 0.8 and 7' 0.4. The left set of data (black) 
corresponds to a trap density p = 0.1 and a slope of 0.378; the right set (red) to 
p — 0.01 and a slope of 0.396. The asymptotic prediction for the slope is 7/2 = 0.4. 
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Table 1. Numerical simulation results for the apparent exponent in the survival 
probability of the particle A. The first five sets of results are for the parameter regime 
7 < 7'. The next four sets are for 7 > 7'. Only the last set is in the regime < 7 < 2/3, 
7' = 1 for which we have no bounding results. In addition we show the x^/d.o.f of the 
fit - d.o.f being the number of degrees of freedom of the fit-, the number of realizations 
performed, Nji, and the logarithm of the minimum value of time used on each fit, 
Mt^l-^J. The statistical error in 9 (one standard deviation) has been computed using 
a jackknife procedure. 



Traps (7) 


Particle(7') p 


9 


7/2 




X^/d.o.f 


Nr 








7 < 7' 










0.4 


0.5 


0.01 


0.146(7) 


0.2 


27.2 


0.58 


460000 






0.1 


0.175(3) 


0.2 


14.47 


0.37 


511733 


0.4 


0.6 


0.01 


0.121(4) 


0.2 


25.12 


0.26 


2050560 






0.1 


0.155(4) 


0.2 


15.74 


1.05 


1398048 


0.4 


0.7 


0.1 


0.093(5) 


0.2 


24.8 


0.81 


6908894 






0.5 


0.190(4) 


0.2 


8.17 


1.04 


2106963 


0.5 


0.9 


0.01 


0.124(2) 


0.25 


18.5 


0.6 


38523723 






0.1 


0.176(4) 


0.25 


11.9 


1.08 


2285460 






0.5 


0.262(6) 


0.25 


6.25 


0.93 


789940 


0.8 


i.U 


0.01 


0.423(2) 


0.4 


12 


1.00 


13626177 






0.1 


0.430(5) 


0.4 


3.55 


0.83 


618793 








7 > 7' 










0.4 


0.4 


0.1 


0.186(2) 


0.2 


15.2 


0.43 


548122 






0.5 


0.210(2) 


0.2 


6.47 


1.05 


425840 


0.8 


0.4 


0.01 


0.396(4) 


0.4 


11.5 


0.4 


111208 






0.1 


0.378(2) 


0.4 


5.58 


0.7 


225796 


0.8 


0.8 


0.01 


0.374(2) 


0.4 


12.14 


1.03 


325275 






0.1 


0.361(3) 


0.4 


8.64 


1.03 


618793 


1.0 


0.5 


0.01 


0.501(2) 


0.5 


7.5 


0.5 


3159970 






< 7 < 2/3, 7' 


= 1 








0.4 


1.0 


0.01 


0.427(6) 


0.2 


10 


1.06 


11952 






0.1 


0.433(7) 


0.2 


5 


1.05 


7740 






0.5 


0.474(7) 


0.2 


3 


1.0 


10000 



the first five sets of results in table [TJ Although (as we will see in the figures) a clear slope 
can be read off the simulation data for ln[— InP(t)], this slope is not independent of p in 
most cases, nor does it yet satisfactorily approach the theoretically predicted asymptotic 
value. We conclude that in this parameter regime asymptotia requires a much longer 
time than we can reasonably simulate, that is, when the particle moves "more easily" 
than the traps, it takes the system a much longer time to behave as it would if the 
particle were simply sitting still. Nevertheless, within this range of parameters there 
are some statements of "quality" than can be made, as shown in figures [5l [6l and [3 
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Figure 4. Simulation results for 7 = 1 and 7' = 0.5. Only one set of data points is 
shown, for p = 0.01. The slope is 0.501. The asymptotic prediction for the slope is 
7/2 = 0.5. 

The most salient point seems to be that asymptotia is reached more readily when both 
7 and 7' are closer to the diffusive case and closer to one another. Thus the results in 
figure [7] are satisfactory (i.e., essentially independent of p and in fair agreement with 
the asymptotic slope) within a reasonable time range. If 7' is close to unity but 7 is too 
small, it is clearly difficult to reach asymptotia, as evidenced in the results of figure El 
If both 7 and 7' are small then even going to extraordinarily long times as shown in 
figure [5] is not sufficient. 

Finally, it is apparent that the simulation results for the "problem case" < 7 < 
2/3, 7' = 1 lead to a slope that seems well defined (see figure [8]) and relatively insensitive 
to p. It is therefore tempting to relate this slope to asymptotic behavior. However, the 
value of the slope does not conffim either of the conjectures put forther earlier, namely, 
that it is perhaps equal to 7/2 or to 1/3. The slope is not particularly close to either 
of these values. On the other hand, the observation does not necessarily disprove the 
conjectures, since it is not clear that asymptotia has been reached. In fact, the observed 
slope is closer to that expected for the short-time behavior of a normally diffusing 
particle surrounded by stationary traps, 9 = 1/2. Thus the question of the asymptotic 
behavior in this regime is still open. 

We have thus found up to this point that our numerical simulations are able 
to conffim the asymptotic prediction for the survival probability of a particle A 
characterized by (sub) diffusion exponent 7' surrounded by traps B characterized by 
exponent 7 provided that 7 < 7', but that it is difficult to do so for 7 > 7'. For the 
regime 0<7<2/3, 7' = ! we have no asymptotic theory, and the numerical results 
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Figure 5. Simulation results for 7 = 0.4 and 7' = 0.5. The left set of data (black in 
colour version) corresponds to a trap density p = 0.1 and a slope of 0.175; the right set 
(red in the colour version) to p = 0.01 and a slope of 0.146. The asymptotic prediction 
for the slope is 7/2 = 0.2. 



3.0 




0.5! I 

2 4 6 8 10 12 14 16 18 20 22 24 26 

ln(t) 

Figure 6. Simulation results for 7 — 0.5 and 7' = 0.9. The left set of data (black) 
corresponds to a trap density p — 0.5 and a slope of 0.262; the middle set (red) to 
p = 0.1 and a slope of 0.176; the right set (blue) to p = 0.01 and a slope of 0.124. The 
asymptotic prediction for the slope is 7/2 = 0.25. 
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Figure 7. Simulation results for 7 = 0.8 and 7' = 1. The left set of data (black) 
corresponds to a trap density p = 0.1 and a slope of 0.430; the right set (red) to 
p = 0.01 and a slope of 0.421. The asymptotic prediction for the slope is 7/2 = 0.4. 



1 1 1 1 1 1 1 1 1 1 1 1 1 1 r 




-2 2 4 6 8 10 12 14 16 

ln(t) 

Figure 8. Simulation results for 7 ~ 0.4 and 7' = 1. The left set of data (black) 
corresponds to a trap density p — 0.5 and a slope of 0.474; the middle set (red) to 
p = 0.1 and a slope of 0.433; the right set (blue) to p = 0.01 and a slope of 0.427. 
There is no asymptotic prediction for the slope in this regime, but two conjectures 
would be that it might be 1/3 or 7/2 — 0.2. 
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Figure 9. Simulation results for the exponential prefactor A in equation ([8]) for 
7 = 7' = 0.4. Broken line: upper bound; Dotted curve: lower bound. Data points are 
for p = 0.1 (red, lower set) and p = 0.5 (black, upper set). 

do not inform us about the validity of conjectured behaviors. 

Having determined the exponent 6 in the survival probability expression ([8]) in some 
parameter regimes, it remains to explore whether we can numerically determine, or at 
least bound, the exponential prefactor A. This turns out to be difficult. In figure [9] we 
show our simulation results for the case 7 = 7' = 0.4 of figure [2l We also show the upper 
bound, which is exactly the asymptotic prediction, cf. compare equations (jl]) and (^^, 
and so appears as a straight (broken) line in the figure. It lies in the lower part of the 
figure because of the minus sign in the ordinate. The lower bound of equation ([5]), which 
only approaches the upper bound asymptotically, is shown as the dotted curve. The 
simulation results are for p = 0.1 and p = 0.5, and fall between the bounds. However, 
we would have to go to times far longer than we are able to in order to ascertain the 
asymptotic prediction. A similar figure, but for only one concentration, p = 0.01, is 
shown in figure [10 for 7 = 1 and 7' = 0.5. Here the upper and lower bounds do 
not even find their rightful relative placements until a time far beyond our simulation 
capabilities, although the simulation results at least point in the right direction. It is, in 
any case, clearly very difficult to determine the prefactor A and even to ascertain that 
it is properly bounded by the theory. 

4. Recap 

In this paper we have made an attempt to assess the validity of the asymptotic 
predictions for the survival probability of a (sub) diffusive particle A characterized by 
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Figure 10. Simulation results for the exponential prefactor A in equation ([5]) for 7 = 1 
and 7' = 0.5. Broken line: upper bound; Dotted curve: lower bound. Data points are 
for p = 0.01. 

exponent 7' surrounded by (sub) diffusive traps of density p characterized by exponent 
7. The prediction is arrived at by obtaining an upper and a lower bound to the survival 
probability that in most parameter regimes converge to one another [19, 20]. This 
asymptotic survival probability in fact turns out to be exactly the upper bound, which 
is calculated under the assumption that particle A remains still. It is thus the case 
that in the parameter regime where this prediction is valid it eventually makes no 
difference whether or not particle A moves; the asymptotic survival probability is entirle 
determined by the motion of the traps. However, when 7' = 1 and the traps are "too 
slow" (0 < 7 < 2/3), the bounds no longer converge even asymptotically, and this 
approach does not lead to a prediction. In other words, it is no longer evident that 
the motion of the particle does not matter. We have proposed two conjectures for this 
regime. One is that in fact the motion of the particle does not matter, as before, but 
our numerical results do not seem to support this assumption. The other relies on 
the fact that for a diffusive particle we know something about the asymptotic survival 
probability at the two extreme points of this interval, namely at 7 = (when the traps 
are stationary) and at 7 = 2/3. In both of these cases the survival probability decays 
as P{t) ~ exp (^—Xt^^^^ (with A known for the former but not for the latter), and so 
one might conjecture a t^^^ dependence in the unknown range. However, this conjecture 
could not be verified either. 

In the regimes where there is an asymptotic prediction, we are able to verify it 
quite clearly when 7 > 7', that is, when the particle moves more slowly or at the same 
pace than the traps. Again, the results indicate that the particle could just as well sit 
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still to reach the same asymptotic survival probability as it does when it moves. Also, 
the "time to asymptotia" is insensitive to the value of 7', but it is shorter when 7 is 
larger and when the density of traps is higher. We also tested our ability to predict 
the asymptotic exponential prefactor A, but find that at best we can show that it lies 
between the correct bounds. At worst, the bounds do not take their rightful places until 
times that we can not reach with our simulations. 

When 7 < 7' the situation is far more difficult, increasingly so with increasing 
difference between the two exponents. It would seem to be necessary to go beyond the 
leading asymptotic term to thoroughly understand the dynamics for these cases. This 
has been done with some measure of success in the purely diffusive problem [40]. 

Our simulation method can not be stretched beyond the times implemented in this 
work. We have been able to answer some questions and ascertain some predictions, 
but not others. To reach the longer times needed to deal with the questions that 
we have not been able answer conclusively will require new simulation optimization 
methods. Such methods have been developed for diffusive particles and traps [37], but 
their generalization to the subdiffusive problem does not appear evident. 
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